function [Z,Zx,Zy] = import_hawaii_dem(july)%filename, domain)
%imports hawaii dem
    %inputs
%filename = name of geotiff dem file, string
%domain = section of dem to import, vector
    %outputs
%Z = ice surface elevation dem, matrix
%Zx,Zy = x,y coordinates of Z
if exist('july','var')
    if ~july
        filename = 'grdn20w156_13_ProjectRaster1.tif';
    else
        filename = 'July2018_Halemaumau_DEM.tif';
    end
else
    filename = 'grdn20w156_13_ProjectRaster1.tif';
end
domain = [-inf,inf,-inf,inf];




[Z, R1] = geotiffread(filename);
Zsize = size(Z);
[Zx,~] = intrinsicToWorld(R1,1:Zsize(2),ones(1,Zsize(2)));
[~,Zy] = intrinsicToWorld(R1,ones(1,Zsize(1)),1:Zsize(1));
Zxi = Zx<=domain(2) & Zx>=domain(1);
Zyi = Zy<=domain(4) & Zy>=domain(3);
Z = Z(Zyi,Zxi);
Zx = Zx(Zxi).';
Zy = Zy(Zyi).';
if Zx(2)<Zx(1)
    Zx = flipud(Zx);
    Z = fliplr(Z);
end
if Zy(2)<Zy(1)
    Zy = flipud(Zy);
    Z = flipud(Z);
end
Z(Z<-9000) = mean(mean(Z(Z>-9000)));
end

